% This version is for chopping experiment. One cycle is background and one
% cycle is data. 
clearvars
close all
%
b.files ={'20230514_100VR_30G_10kHz_REMPI_scan_3'};
%{'20230122_100VR_30G_10kHz_scan_4'};
%{'20221229_100VR_30G_Rydberg_scan_17','20221228_100VR_30G_Rydberg_scan_7','20221229_100VR_30G_Rydberg_scan_16'}};{'20230109_100VR_30G_Rydberg_scan_1'};{{'20230105_100VR_30G_Rydberg_scan_9','20230105_100VR_30G_Rydberg_scan_7'}};{{'20221025_2000V_2G_Rydberg_scan_2','20221025_2000V_2G_Rydberg_scan_3','20221025_2000V_2G_Rydberg_scan_4','20221025_2000V_2G_Rydberg_scan_5','20221026_2000V_2G_Rydberg_scan_2','20221026_2000V_2G_Rydberg_scan_3','20221026_2000V_2G_Rydberg_scan_5','20221026_2000V_2G_Rydberg_scan_6'}};{'20221020_2000V_2G_Rydberg_scan_14'};{{'20220910_2000V_2G_Rydberg_scan_2','20220910_2000V_2G_Rydberg_scan_1'}};{{'20220908_2000V_2G_Rydberg_scan_19','20220908_2000V_2G_Rydberg_scan_18'}};{'20220831_2000V_2G_Rydberg_scan_7','20220901_2000V_2G_Rydberg_scan_1'};{{'20220627_400VR_0G_Rydberg_scan_2','20220628_400VR_0G_Rydberg_scan_1','20220702_400VR_0G_Rydberg_scan_1'}};{'20220722_400VR_0G_Rydberg_scan_2'};{'20220723_400VR_0G_Rydberg_scan_1'};{'20220709_400VR_0G_Rydberg_scan_1','20220712_400VR_0G_Rydberg_scan_1', '20220709_400VR_0G_Rydberg_scan_1','20220714_400VR_0G_Rydberg_scan_1','20220714_400VR_0G_Rydberg_scan_2','20220714_400VR_0G_Rydberg_scan_3'};{{'20220712_400VR_0G_Rydberg_scan_1', '20220709_400VR_0G_Rydberg_scan_1'}};{'20220708_400VR_0G_Rydberg_scan_6'};{'20220707_400VR_0G_Rydberg_scan_9',{'20220707_400VR_0G_Rydberg_scan_9','20220707_400VR_0G_Rydberg_scan_9'}};

b.TOF = {[48930, 69196, 72150, 87180, 99971, 102045, 113167, 122000]*1.082 - 19500  ... % 10kHz
     [48930, 69196, 72150, 87180, 99971, 102045, 113167, 122000]*1.082 - 19500   % 10kHz
    };
b.TOFWindow = 1800;1200; %[ns]
b.species = [0, 0,1, 1, 0,0,0,0];
b.cyclesCut = {[],[],[]}; {[], [50 57;82 84]}; {[],[],[]};{[0 29], [50 57;82 84],[2 7;74 79;88 93],[22 27;28 31;44 51;54 67],[12 19;41 51],[],[],[26 31],[72 75]};                {[8 16;76 78 ;90 97; 105 109],[],[26 31],[72 75]}; {[0 0],[],[],[0,0], [0,0]};{[100 215],[0 0],[0,0],[0,0],[0,0],[0,0]};{[260,265;464,1429],[0,0]};{[700,981],[0,0],[0,0],[0,0]}; % cycles to trash
% b.cyclesCut =  {[260,265;464,1429],[0,0]};

b.maxShot = {2*10000,2*10000,2*5000,2*5000,2*5000,2*5000,2*5000,2*5000,2*5000,2*5000};
b.shotCut = {[000,10000],[000,10000], [0,10000],[1000,5000],[1000,5000],[1000,5000], [1000,5000],[1000,5000],[1000,5000],[1000,5000],}; % shots to keep
b.xCut = {[],[],[],[],[],[],[],[]}; % x variables to keep 
% b.xCut = {[0 1000],[0 1000],[0 1000],[0 1000],[0 1000],[0 1000],[460 490],[460 490]}; % x variables to keep 
b.xShift = {[0],[0],[0],[0],[0],[0],[0],[0]};
b.xShift = {[0],[0],[0],[0],[0],[0],[0],[0]};
b.normalize = 0;
b.rewriteXvar = 0;
b.speciesNames = {'K','K_2','Rb','KRb','K_2Rb','Rb_2','KRb_2','K_2Rb_2'};
b.xlabel = {'420 frequency 713288.510GHz + x MHz'};{'670 frequency 449846 GHz + x [MHz]'}; {'670 delay 43.02 us + x ns'};{'ARP duration (ms)'};{'Detection delay (ms)'}; {'670 frequency 449850 GHz + x (MHz)'};{'11 to 22 uwave duration (pi pulse is 37.1 us) (us)'};   {'670 delay 43.02 us + x [ns]'};{'pi pulse duration us'};{'420 42.46us + x ns, 670 43.03us + x ns'}; {'670: 449844.24 GHz + x [MHz]'};{'Delay time between 1,1 to 2,2 pi pulse and ionization (ms)'};{'670 frequency (MHz) + 449850 GHz'};{'UV delay (us)'};{'670 frequency (MHz) + 449850 GHz'}; {'Wait time after pi pulse (us)'};{'420 frequency 713288 GHz + x [MHz]'};;{'420 frequency (MHz) + 713288 GHz'};
b.comments = {'|1,-1> (ARP at 30 G) ','|1 1> ','|2,2> (pi pulse @ 2G) + |0,0,-4,1/2>','|1,1> + |0,0,-4,1/2> 2G','|1,1> + |0,0,-4,3/2>'};{'|2,2> (ARP at 542G) + N=0', '|1,1> + N=0',' |2,2> (pi pulse at 2G) + N=0',}; {'with |2,2> atoms','|2,2> atoms optically killed','19k Rb atoms in |2,2>'};{'Rb number 12k','Rb number 3k'};{'IR resonant','IR off-resonant'};{'35D Rb + Rb','35D Rb + N=0 KRb','35D Rb + |1,0,-4,1/2> KRb','','',''};{'Rb-Rb 20220704','Rb-Rb 20220705'};{'new data','Rb-KRb combined data'};
b.chop.flag = 1; % 0: no chop; 1: background first then data; 2: data first and then data
b.chop.T = 2;
b.var = 1; % use shot to shot fluctuation as a representation of errorbar


%% load all the data sets
DataSets = {};
set_num = length(b.files);
for i = 1:set_num    
    fi = b.files{i}; 
    if iscell(fi) && length(fi) > 1
        for j = 1:length(fi)
            p1 = load_data(fi{j}); 
            DataSets{end+1} = p1;
            clear p1
        end
    else
        p = load_data(fi);
        DataSets{end+1} = p;
    end    
end

%% process using modualized functions
% processed data stored in a cell array, each element is a matrix
dataAll = {}; % all filtered data
pdAll = {}; % all processed data, 2d cell array, the first dimension is for groups in chopping, second dimension is for different datasets

for i = 1:length(DataSets)
    [filtered_data] =  filter_by_shot_and_cycle(DataSets,b,i);
    dataAll{i} = filtered_data;
    
    for j = 1:b.chop.T
        [pd] = clean_up_by_xVar_and_species(filtered_data{j},b,i);
        pdAll{j}{i} = pd;
    end
end

%% combine data sets
cpdAll = {}; % combined processed datasets {pd1,pd2,pd3}; pd1 dim: (length(xVar),8); 8 species
for j = 1:b.chop.T
    cpdAll{end+1} = combine_data(pdAll{j},b);
end
%% plotting
% plot TOF histogram
% scrsz = get(groot, 'Screensize');
% h5 = figure();
% set(h5, 'Position', [scrsz(3) 1 scrsz(3) scrsz(4)])
% set(h5,'Color','w')
% for i = 3:3
%     plot_TOF(dataAll{i},b,i,h5)
% end

% plot differnet species 
% plot_by_x_var(cpdAll{1},b)
% plot_by_x_var(cpdAll{1},b)
plot_background(cpdAll,b)



% % plot different species on the same plot
% plot_overlap_species(cpdAll,b,1)

% % plot combined positions
% cPos = plot_combined_position(dataAll,b);


%% functions
function p = load_data(filename)
    cd('\\krb1-pc\d\Data analysis scripts\MCP Image Analysis\MCP Shot by Shot\lmf2txt\');
    folderpath = '\\KRB1-PC\d\Data analysis scripts\MCP Image Analysis\MCP Shot by Shot\lmf2txt\';

    filepath = [folderpath, filename];
    if (~exist([filepath,'.dat'],'file'))
        error([filepath,'.dat',' needs to be created using mainMCPtxtAnalysisTiming.m first!'])
    else
        p = dlmread([filepath,'.dat']);
    end
end

function [filtered_data] = filter_by_shot_and_cycle(DataSets,b,ind)
    % also do rotation here
    realHit = DataSets{ind};
    cyclesCut = b.cyclesCut{ind};
    flag_rewrite_xVar = b.rewriteXvar;
    % Name columns of realHit
    TOFAllShots = realHit(:,1);
    xPosAllShots = realHit(:,2);
    yPosAllShots = realHit(:,3);
    shotNumAllShots = realHit(:,4);
    xVarAllShots1 = realHit(:,5);
    xVarAllShots1014 = realHit(:,8);
    
    % Calculate the total number of experimental cycles
    totCycleNum = 1;
    cycleNumAllShots = zeros(length(realHit),1);
    cycleNumAllShots(1) = 1;

    % assign a cycle number to every shot. 
    for i = 2:length(shotNumAllShots)
        if shotNumAllShots(i) < shotNumAllShots(i-1)
            totCycleNum = totCycleNum + 1;
        end
        cycleNumAllShots(i) = totCycleNum;
    end
    
    [numCutPair,~] = size(cyclesCut);
    % Comments about "cyclesCut"
    % Each pair of values specify an interval of cycles to be cut out of the data analysis
    % Example: to cut from shots 100 - 200 and 300 - 400, use the format cyclesCut = [100 200;300 400];
    % To not cut anything, put [0 0]
    cycleNumInd = 1:length(cycleNumAllShots);
    for i = 1:numCutPair
        tmp = find((cycleNumAllShots <= cyclesCut(i,1))|(cycleNumAllShots >= cyclesCut(i,2))); % find indices of shots that will survive the cyclesCut
        cycleNumInd = intersect(cycleNumInd,tmp);
    end
    
    maxShot = b.maxShot{ind};
    shotNumStart = b.shotCut{ind}(1); % 300
    shotNumEnd =  b.shotCut{ind}(2); %maxShot./2
      
    % define data vector filtered by shot number and cycle number
    shotNumIndAll = find((shotNumAllShots >= shotNumStart)&(shotNumAllShots <= shotNumEnd)); % returns indices   
    filterNumIndAll = intersect(shotNumIndAll, cycleNumInd); % returns the common element of shotNumIndAll and cycleNumInd
    
    % group data (for chopping experiment)
    filtered_data = {};
    T = b.chop.T;
    for i = 1:T
        if b.chop.flag ~= 0
%             group_ind = find(mod(cycleNumAllShots,T)==T-i); % This is
%             wrong 
            if i==T
                group_ind = find(mod(cycleNumAllShots,T)==0);
            else
                group_ind = find(mod(cycleNumAllShots,T)==i);
            end
            filterNumIndAll_temp = intersect(group_ind,filterNumIndAll);
            
        else
            filterNumIndAll_temp = filterNumIndAll;
        end
        data.TOF = TOFAllShots(filterNumIndAll_temp);
        xPosAll = xPosAllShots(filterNumIndAll_temp);
        yPosAll = yPosAllShots(filterNumIndAll_temp);
        data.shotNum = shotNumAllShots(filterNumIndAll_temp);
        data.cycleNum = cycleNumAllShots(filterNumIndAll_temp);
        data.xVar = xVarAllShots1(filterNumIndAll_temp);
        if flag_rewrite_xVar
            data.xVar = xVarAllShots1014(filterNumIndAll_temp);
        end
        data.totCycleNumEff = length(unique(data.cycleNum));
        data.totCycleNum = length(unique(cycleNumAllShots));
        
        % match the number of experiments and background
        shotNumBkgdStart = maxShot./2 + 1;
        shotNumBkgdEnd = shotNumBkgdStart + shotNumEnd - shotNumStart ;
        
        % Define bkgd data vectors filtered by shot number and cycle number
        shotNumIndBkgd = find((shotNumAllShots >= shotNumBkgdStart)&(shotNumAllShots <= shotNumBkgdEnd));
        filterNumIndBkgd = intersect(shotNumIndBkgd, cycleNumInd);
        
        bkgd.TOF = TOFAllShots(filterNumIndBkgd);
        xPosBkgd = xPosAllShots(filterNumIndBkgd);
        yPosBkgd = yPosAllShots(filterNumIndBkgd);
        bkgd.shotNum = shotNumAllShots(filterNumIndBkgd);
        bkgd.cycleNum = cycleNumAllShots(filterNumIndBkgd);
        bkgd.xVar = xVarAllShots1(filterNumIndBkgd);
        if flag_rewrite_xVar
            bkgd.xVar = xVarAllShots1014(filterNumIndBkgd);
        end
        
        % Rotation
        theta = 45;
        data.xPos = cos(theta/180*pi)*xPosAll - sin(theta/180*pi)*yPosAll;
        data.yPos = sin(theta/180*pi)*xPosAll + cos(theta/180*pi)*yPosAll;
        bkgd.xPos = cos(theta/180*pi)*xPosBkgd - sin(theta/180*pi)*yPosBkgd;
        bkgd.yPos = sin(theta/180*pi)*xPosBkgd + cos(theta/180*pi)*yPosBkgd;
        
        % get TOF histogram
        speciesMasses = [40,80,87,127,167,174,214,254];
        TOFbinSizeAll = 50; %in ns
        TOFWindow = b.TOFWindow; %1200
        TOFBinSizeSpecies = [1, 10, 1, 1, 1, 10, 1, 1];
        
        TOFMin = 0;
        TOFMax = 100000; % [ns] was 100000 ns
        TOFEdgesAll = TOFMin:TOFbinSizeAll:TOFMax;
        
        data.tofCounts = histcounts(data.TOF,TOFEdgesAll)';
        bkgd.tofCounts = histcounts(bkgd.TOF,TOFEdgesAll)';
        data.t = (TOFEdgesAll(1:length(TOFEdgesAll)-1) + TOFEdgesAll(2:length(TOFEdgesAll)))./2;
        data.t = data.t';
        bkgd.t = data.t;
        
        fd.data = data;
        fd.bkgd = bkgd;
        filtered_data{end+1} = fd;
    end
end

function [pd] = clean_up_by_xVar_and_species(filtered_data,b,ind)
    
    data = filtered_data.data;
    bkgd = filtered_data.bkgd;

    MCPRadius = 40;
    xMin = -MCPRadius;
    xMax = MCPRadius;
    yMin = -MCPRadius;
    yMax = MCPRadius;
    
    xCentROISpecies = [12.95, 5.024, -10, 1.37, 0, 10, -2, -1.5]; %[12.95, 5.024, 5.05, 1.37, 0, -0.72, -0.5, -1.5];
    yCentROISpecies = [-0.85, -0.5887, -0.85, -0.85, 0, -0.85, -0.85, -0.85];
    windowROISpecies = [80, 80, 80, 80, 80, 80, 80, 80];
    
    
    
    theta = 45;
    xCentROISpecies = [-16, 5.85, -10, 2.25, 0, 0.65, 0, 0];
    yCentROISpecies = [-1, -0.85, -1, -0.85, 0, -0.85, 0, 0];
    windowROISpecies = [10, 10, 80, 80, 10, 10, 80, 10];
    % Calculatex ROI windows for different species
    xMinROISpecies = xCentROISpecies - windowROISpecies./2;
    xMaxROISpecies = xCentROISpecies + windowROISpecies./2;
    yMinROISpecies = yCentROISpecies - windowROISpecies./2;
    yMaxROISpecies = yCentROISpecies + windowROISpecies./2;
    
    xVarList1 = unique(data.xVar);
    totCountsSpeciesxVar = zeros(length(xVarList1),length(b.speciesNames));
    totCountsSpeciesBkgdxVar = zeros(length(xVarList1),length(b.speciesNames));
    cycleNumxVar = zeros(size(xVarList1));

    TOFWindow = b.TOFWindow; %1200
    TOFAll = data.TOF;
    TOFs = b.TOF{ind};
    TOFBkgd = bkgd.TOF;
    for j = 1:length(xVarList1)
        for i = 1:length(b.speciesNames)
            TOFSpeciesInd = find((TOFAll >= TOFs(i) - TOFWindow/2) & (TOFAll <= TOFs(i) + TOFWindow/2) ...
                & (data.xPos >= xMinROISpecies(i)) & (data.xPos <= xMaxROISpecies(i)) ...
                & (data.yPos >= yMinROISpecies(i)) & (data.yPos <= yMaxROISpecies(i)) ...
                & (data.xVar == xVarList1(j)));
            TOFSpeciesBkgdInd = find((TOFBkgd >= TOFs(i) - TOFWindow/2) & (TOFBkgd <= TOFs(i) + TOFWindow/2) ...
                & (bkgd.xPos >= xMinROISpecies(i)) & (bkgd.xPos <= xMaxROISpecies(i)) ...
                & (bkgd.yPos >= yMinROISpecies(i)) & (bkgd.yPos <= yMaxROISpecies(i)) ...
                & (bkgd.xVar == xVarList1(j)));

            TOFSpecies = TOFAll(TOFSpeciesInd);
            totCountsSpeciesxVar(j,i) = length(TOFSpecies);

            TOFSpeciesBkgd = TOFBkgd(TOFSpeciesBkgdInd);
            totCountsSpeciesBkgdxVar(j,i) = length(TOFSpeciesBkgd);
            
            if b.var % use shot to shot fluctuation to assign errorbar 
            end

        end
        xVarListInd = find(data.xVar == xVarList1(j));
        cycleNumxVar(j) = length(unique(data.cycleNum(xVarListInd)));
    end


    totCountsSpeciesRealxVar = totCountsSpeciesxVar - totCountsSpeciesBkgdxVar;
    totCountsSpeciesRealxVar_Err = sqrt(totCountsSpeciesxVar + totCountsSpeciesBkgdxVar);
    totCountsSpeciesRealxVarNorm = totCountsSpeciesRealxVar./cycleNumxVar;
    %     totCountsSpeciesRealxVarNorm_Err = totCountsSpeciesRealxVar_Err./totCountsSpeciesRealxVar.*totCountsSpeciesRealxVarNorm;
    totCountsSpeciesRealxVarNorm_Err = totCountsSpeciesRealxVar_Err./cycleNumxVar;
    pd.xVar = xVarList1;
    pd.data = totCountsSpeciesRealxVarNorm;
    pd.err = totCountsSpeciesRealxVarNorm_Err;
    pd.rep = cycleNumxVar;

    if ~isempty(b.xCut{ind})
        g = find(xVarList1>=b.xCut{ind}(1) & xVarList1 <= b.xCut{ind}(2));
        pd.xVar = xVarList1(g);
        pd.data = totCountsSpeciesRealxVarNorm(g,:);
        pd.err = totCountsSpeciesRealxVarNorm_Err(g,:);
        pd.rep = cycleNumxVar(g,:);
    end

end

function  [cPos] = plot_combined_position(dataAll,b)
    cpx = [];
    cpy = [];
    cTOF = [];
    for i = 1: length(dataAll)
        cpx = cat(1,cpx, dataAll{i}.data.xPos);
        cpy = cat(1, cpy, dataAll{i}.data.yPos);
        cTOF = cat(1, cTOF, dataAll{i}.data.TOF);
    end
    figure(2)
    h2 = figure(); % Spatial figure
    scrsz = get( groot, 'Screensize' );
    set(h2, 'Position', [scrsz(3) 1 scrsz(3) scrsz(4)])
    %     set(h2, 'Position', [scrsz(3) 1 scrsz(3) scrsz(4)])
    %set colormap for MCP spatial image
    cmap = jet(5000);
    cmap(1,:) = zeros(1,3);
    set(gcf,'Color','w')
    MCPRadius = 40;
    xMin = -MCPRadius;
    xMax = MCPRadius;
    yMin = -MCPRadius;
    yMax = MCPRadius;
    theta = 45;
    xCentROISpecies = [12.95, 5.024, 10, 1.37, 0, 10, -2, -1.5]; %[12.95, 5.024, 5.05, 1.37, 0, -0.72, -0.5, -1.5];
    yCentROISpecies = [-0.85, -0.5887, -0.85, -0.85, 0, -0.85, -0.85, -0.85];
    windowROISpecies = [80, 10, 80, 80, 10, 80, 80, 10];
    spatialBinSizeSpecies = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1];
    %% Calculated ROI windows for different species
    xMinROISpecies = xCentROISpecies - windowROISpecies./2;
    xMaxROISpecies = xCentROISpecies + windowROISpecies./2;
    yMinROISpecies = yCentROISpecies - windowROISpecies./2;
    yMaxROISpecies = yCentROISpecies + windowROISpecies./2;
    
    TOFs = b.TOF{1};
    TOFWindow = b.TOFWindow; 
    cPos = {};
    for i = 1:length(b.speciesNames)
        TOFSpeciesInd = find((cTOF >= TOFs(i) - TOFWindow/2) & (cTOF <= TOFs(i) + TOFWindow/2) ...
            & (cpx >= xMinROISpecies(i)) & (cpx <= xMaxROISpecies(i)) ...
            & (cpy >= yMinROISpecies(i)) & (cpy <= yMaxROISpecies(i)));
        xPosSpecies = cpx(TOFSpeciesInd);
        yPosSpecies = cpy(TOFSpeciesInd);
        xEdgesSpecies = xMin:spatialBinSizeSpecies(i):xMax;
        yEdgesSpecies = yMin:spatialBinSizeSpecies(i):yMax;
        xSpecies = (xEdgesSpecies(1:length(xEdgesSpecies)-1) + xEdgesSpecies(2:length(xEdgesSpecies)))/2;
        ySpecies = (yEdgesSpecies(1:length(yEdgesSpecies)-1) + yEdgesSpecies(2:length(yEdgesSpecies)))/2;
        
        xPosSpeciesRot = xPosSpecies;
        yPosSpeciesRot = yPosSpecies;
        spatialCountsSpeciesRot = transpose(histcounts2(xPosSpeciesRot, yPosSpeciesRot, xEdgesSpecies, yEdgesSpecies));
        %     spatialCountsSpeciesRot = imrotate(spatialCountsSpecies,theta,'crop');
        xCentSpecies(i) = mean(xPosSpeciesRot);
        yCentSpecies(i) = mean(yPosSpeciesRot);
        if b.species(i)
            if length(nonzeros(b.species)) <= 4
                subplot(1,length(nonzeros(b.species)),find(find(b.species == 1)==i))
            else
                subplot(2,ceil(length(nonzeros(b.species))/2),find(find(b.species == 1)==i))
            end
            colormap(cmap);
            
            if max(max(spatialCountsSpeciesRot)) < 2
                spatialCountsSpeciesRot(1,1) = 2;
            end
            
            %         spatialCountsSpeciesRot(spatialCountsSpeciesRot>20)=20;
            imagesc(xSpecies, ySpecies, spatialCountsSpeciesRot);
            
            cp.x = xSpecies;
            cp.y = ySpecies;
            cp.z = spatialCountsSpeciesRot;
            cp.xlist = xPosSpeciesRot;
            cp.ylist = yPosSpeciesRot;
            cPos{end+1} = cp;
            %         imagesc(spatialCountsSpeciesRot);
            set(gca,'YDir','normal','XDir','normal')
            hold on
            k = linspace(0,2*pi);
            plot(MCPRadius.*cos(k),MCPRadius.*sin(k),'w');
            
            plot([xMinROISpecies(i),xMinROISpecies(i)],[yMinROISpecies(i),yMaxROISpecies(i)],'w')
            hold on
            plot([xMaxROISpecies(i),xMaxROISpecies(i)],[yMinROISpecies(i),yMaxROISpecies(i)],'w')
            hold on
            plot([xMinROISpecies(i),xMaxROISpecies(i)],[yMinROISpecies(i),yMinROISpecies(i)],'w')
            hold on
            plot([xMinROISpecies(i),xMaxROISpecies(i)],[yMaxROISpecies(i),yMaxROISpecies(i)],'w')
            hold on         
            
            axis square
            hold off
            colorbar('east','Color','w')
            xlabel('x (mm)', 'FontSize', 14)
            ylabel('y (mm)', 'FontSize', 14)
            
            xlim([xMinROISpecies(i), xMaxROISpecies(i)])
            ylim([yMinROISpecies(i), yMaxROISpecies(i)])
            text(xMinROISpecies(i) + 0.01*(xMaxROISpecies(i) - xMinROISpecies(i)), yMinROISpecies(i) + 0.15*(yMaxROISpecies(i) - yMinROISpecies(i)), ['N_{ROI} = ',num2str(length(TOFSpeciesInd))], 'FontSize', 14, 'Color','w')
            text(xMinROISpecies(i) + 0.01*(xMaxROISpecies(i) - xMinROISpecies(i)), yMinROISpecies(i) + 0.95*(yMaxROISpecies(i) - yMinROISpecies(i)), [b.speciesNames{i},'^+'], 'FontSize', 14, 'Color','w')
            
            %             title(['MCP Image of shot # ',num2str(shotNumStart),'-',num2str(shotNumEnd)], 'Interpreter', 'none')
            
        end
    end
    
end

function [cpd] = combine_data(pd,b)
	cpd = {};
    M = {};
    k = 1;
    for i = 1:length(b.files)
        fi = b.files{i};
        if iscell(fi) 
            c = [];
            for j = 1:length(fi)
                c(end+1) = k;
                k = k+1;
            end
            M{end+1} = c;
        else
            M{end+1} = k;
            k = k+1;
        end
    end    
    for i = 1:length(M)      
        if length(M{i})>1
            xVar = [];
            data = [];
            err = [];
            rep = [];
            for j=1:length(M{i})
                ind = M{i}(j);
                xVar = cat(1,xVar,pd{ind}.xVar);
                data = cat(1,data,pd{ind}.data);
                err = cat(1,err,pd{ind}.err);
                rep = cat(1,rep,pd{ind}.rep);
                
                % to calculate the error bar of the combined data
                % first undo the normalization of data by cycle number
                data_tot = data.*rep;
                err_tot = err.*rep;
            end
            
            % findgroups and splitapply
            G = findgroups(xVar);
            xVar = splitapply(@mean, xVar, G);
            rep = splitapply(@sum, rep, G);
            data_tot = splitapply(@(x)sum(x,1), data_tot, G);
            err_tot = splitapply(@(x) sqrt(sum(x.^2,1)), err_tot, G);
            data = data_tot./rep;
            err = err_tot./rep;
            
            cpd1.xVar = xVar;
            cpd1.data = data;
            cpd1.err = err;
            cpd1.rep = sum(rep);
            cpd{i} = cpd1;
        else
            ind = M{i};
            cpd{i} = pd{ind};
            cpd{i}.rep = sum(cpd{i}.rep);
        end
    end  
    
    if b.normalize
        for i = 1:length(cpd)
            pdt = cpd{i};
            for j=1:8
                coeff = max(pdt.data(:,j));
                pdt.data(:,j) = pdt.data(:,j)/coeff;
                pdt.err(:,j) = pdt.err(:,j)/coeff;
            end
            cpd{i} = pdt;
        end
    end
end

function plot_TOF(filtered_data,b,ind,h5)
    data = filtered_data.data;
    bkgd = filtered_data.bkgd;
    figure(h5);
    hold on
    subplot(3,1,1)
    TOFTimeAll = data.t;
    TOFCountsAll = data.tofCounts;
    TOFCountsBkgd = bkgd.tofCounts;
    maxShot = b.maxShot{ind};
    shotNumStart = b.shotCut{ind}(1); % 300
    shotNumEnd =  b.shotCut{ind}(2); %maxShot./2
    totCycleNum = data.totCycleNum;
    totCycleNumEff = data.totCycleNumEff;
    TOFs = b.TOF{ind}; TOFWindow = b.TOFWindow;
    shotNumBkgdStart = maxShot./2 + 1;
    shotNumBkgdEnd = shotNumBkgdStart + shotNumEnd - shotNumStart ;
   
    colorCodes = {'k','r','k','k','g','b','m',[30 144 255]./255};
    
    plot(TOFTimeAll,TOFCountsAll,'k')
%     hold off
    TOFMin = 0;
    TOFMax = 100000;
    xlim([TOFMin,TOFMax])
    ylim([0,1.2*max(TOFCountsAll)])
    grid on
    grid minor
    set(gca,'FontSize',14);
    ylabel('Ion counts', 'FontSize', 16)
    title([b.files{ind}, ', TOF spectrum of shot # ',num2str(shotNumStart),'-',num2str(shotNumEnd)], 'Interpreter', 'none')
    text(TOFMin, 1.1*max(TOFCountsAll), [num2str(totCycleNum), ' total cycles'], 'FontSize', 16);
    text(TOFMin, 0.95*max(TOFCountsAll), [num2str(totCycleNumEff), ' useful cycles'], 'FontSize', 16);
    text(TOFMin, 0.80*max(TOFCountsAll), [num2str(maxShot), ' shots/cycle'], 'FontSize', 16);
    
    for i = 1:length(b.speciesNames)
        if b.species(i)
            patch([TOFs(i)-TOFWindow/2 TOFs(i)-TOFWindow/2 TOFs(i)+TOFWindow/2 TOFs(i)+TOFWindow/2], ...
                [-1 1.3*max(TOFCountsAll) 1.3*max(TOFCountsAll) -1], ...
                colorCodes{i}, 'EdgeColor', [1 1 1])
            alpha 0.20
%             if isequal(b.speciesNames{i},'K') || isequal(b.speciesNames{i},'K_2') || isequal(b.speciesNames{i},'K_2Rb')
%                 text(TOFs(i)+TOFWindow/2, 1.1*max(TOFCountsAll), [b.speciesNames{i},'^+','(N = ',num2str(totCountsSpecies(i)),')'], 'FontSize', 14, 'Color', colorCodes{i})
%             else
%                 text(TOFs(i)+TOFWindow/2, 1.0*max(TOFCountsAll), [b.speciesNames{i},'^+','(N = ',num2str(totCountsSpecies(i)),')'], 'FontSize', 14, 'Color', colorCodes{i})
%             end
            if isequal(b.speciesNames{i},'K') || isequal(b.speciesNames{i},'K_2') || isequal(b.speciesNames{i},'K_2Rb')
                text(TOFs(i)+TOFWindow/2, 1.1*max(TOFCountsAll), [b.speciesNames{i},'^+'], 'FontSize', 14, 'Color', colorCodes{i})
            else
                text(TOFs(i)+TOFWindow/2, 1.0*max(TOFCountsAll), [b.speciesNames{i},'^+'], 'FontSize', 14, 'Color', colorCodes{i})
            end
        end
    end
    
    subplot(3,1,2)
    plot(TOFTimeAll,TOFCountsAll,'k')
%     hold off
    xlim([TOFMin,TOFMax])
    ymaxzoom = 100;
%     ymaxzoom = 1.2*max([totCountsK_2, totCountsRb_2, totCountsK_2Rb,totCountsKRb_2,totCountsK_2Rb_2]); % Zoom in over the y-axis (ion counts) of the TOF spectrum
%     ymaxzoom = 1.2*max([totCountsK_2, totCountsRb_2]);
    ylim([0 ymaxzoom])
    grid on
    grid minor
    set(gca,'FontSize',14);
    % xlabel('Time-of-flight (ns)', 'FontSize', 16)
    ylabel('Ion counts', 'FontSize', 16)
    title([b.files{ind}, ', TOF spectrum of shot # ',num2str(shotNumStart),'-',num2str(shotNumEnd)], 'Interpreter', 'none')
%     text(TOFMin, 1.1*max(TOFCountsAll), [num2str(totCycleNum), ' cycles'], 'FontSize', 16);
%     text(TOFMin, 0.95*max(TOFCountsAll), [num2str(maxShot), ' shots/cycle'], 'FontSize', 16);
    
    for i = 1:length(b.speciesNames)
        if b.species(i)
            patch([TOFs(i)-TOFWindow/2 TOFs(i)-TOFWindow/2 TOFs(i)+TOFWindow/2 TOFs(i)+TOFWindow/2], ...
                [-1 1.3*max(TOFCountsAll) 1.3*max(TOFCountsAll) -1], ...
                colorCodes{i}, 'EdgeColor', [1 1 1])
            alpha 0.20
%             if isequal(b.speciesNames{i},'K') || isequal(b.speciesNames{i},'K_2') || isequal(b.speciesNames{i},'K_2Rb')
%                 text(TOFs(i)+TOFWindow/2, 0.95*ymaxzoom, [b.speciesNames{i},'^+','(N = ',num2str(totCountsSpecies(i)),')'], 'FontSize', 14, 'Color', colorCodes{i})
%             else
%                 text(TOFs(i)+TOFWindow/2, 0.85*ymaxzoom, [b.speciesNames{i},'^+','(N = ',num2str(totCountsSpecies(i)),')'], 'FontSize', 14, 'Color', colorCodes{i})
%             end
            if isequal(b.speciesNames{i},'K') || isequal(b.speciesNames{i},'K_2') || isequal(b.speciesNames{i},'K_2Rb')
                text(TOFs(i)+TOFWindow/2, 0.95*ymaxzoom, [b.speciesNames{i},'^+'], 'FontSize', 14, 'Color', colorCodes{i})
            else
                text(TOFs(i)+TOFWindow/2, 0.85*ymaxzoom, [b.speciesNames{i},'^+'], 'FontSize', 14, 'Color', colorCodes{i})
            end
        end
    end
    
    subplot(3,1,3)
    plot(TOFTimeAll,TOFCountsBkgd,'k')
%     hold off
    xlim([TOFMin,TOFMax])
    ylim([0 ymaxzoom])
    grid on
    grid minor
    set(gca,'FontSize',14);
    xlabel('Time-of-flight (ns)', 'FontSize', 16)
    ylabel('Ion counts', 'FontSize', 16)
    title([b.files{ind}, ', TOF spectrum of shot # ',num2str(shotNumBkgdStart),'-',num2str(shotNumBkgdEnd)], 'Interpreter', 'none')
    
    
    for i = 1:length(b.speciesNames)
        if b.species(i)
            patch([TOFs(i)-TOFWindow/2 TOFs(i)-TOFWindow/2 TOFs(i)+TOFWindow/2 TOFs(i)+TOFWindow/2], ...
                [-1 1.3*max(TOFCountsAll) 1.3*max(TOFCountsAll) -1], ...
                colorCodes{i}, 'EdgeColor', [1 1 1])
            alpha 0.20
%             if isequal(b.speciesNames{i},'K') || isequal(b.speciesNames{i},'K_2') || isequal(b.speciesNames{i},'K_2Rb')
%                 text(TOFs(i)+TOFWindow/2, 0.95*ymaxzoom, [b.speciesNames{i},'^+','(N = ',num2str(totCountsSpeciesBkgd(i)),')'], 'FontSize', 14, 'Color', colorCodes{i})
%             else
%                 text(TOFs(i)+TOFWindow/2, 0.85*ymaxzoom, [b.speciesNames{i},'^+','(N = ',num2str(totCountsSpeciesBkgd(i)),')'], 'FontSize', 14, 'Color', colorCodes{i})
%             end
            if isequal(b.speciesNames{i},'K') || isequal(b.speciesNames{i},'K_2') || isequal(b.speciesNames{i},'K_2Rb')
                text(TOFs(i)+TOFWindow/2, 0.95*ymaxzoom, [b.speciesNames{i},'^+'], 'FontSize', 14, 'Color', colorCodes{i})
            else
                text(TOFs(i)+TOFWindow/2, 0.85*ymaxzoom, [b.speciesNames{i},'^+'], 'FontSize', 14, 'Color', colorCodes{i})
            end
        end
    end
end

function plot_by_x_var(cpd,b)
    h7 = figure(); 
    scrsz = get(groot, 'Screensize');
    set(h7, 'Position', [scrsz(3) + scrsz(3)/2 1 scrsz(3)/2 scrsz(4)])
    set(h7,'Color','w')
    species = b.species;
    pinks = {'#860454','#c51b7c','#dc75ab','#f0b7da','#FF00FF','#ff81c0'};
    greens = {'#276418','#4d921e','#7fbc42','#b9e084','#15b01a','#77AC30'};
    purples = {'#410149','#762a84','#9b6fac','#c1a5cd','#800080','#7e1e9c'};
    browns = {'#523107','#7f3c0a','#bf812c','#e2c17e','#a0522d','#A2142F'};
    pools = {'#003d2e','#01665e','#379692','#81cdc1','#13eac9','#04d8b2'};
    reds = {'#68011d','#b5172f','#d75f4e','#f7a580','#ff0000','#e50000'};
    blues = {'#062e61','#2265ad','#4295c1','#93c5dc','#40e0d0','#00FFFF'};
    b.species = [0, 0, 1, 1, 0, 1 , 1, 0];
    colorCodes = {reds,pools,purples,greens,pinks,blues,browns,pinks};
    for i = 1:length(b.speciesNames)
        if species(i)
            figure(h7)
            hold on
            if length(nonzeros(species)) <= 4
                subplot(length(nonzeros(species)),1,find(find(species == 1)==i))
                hold on
            else
                subplot(ceil(length(nonzeros(species))/2),2, find(find(species == 1)==i))
                hold on
            end
            
            ymin = [];
            ymax = [];
            for k = 1:length(cpd)
                xShift = b.xShift{k};
                xlist = cpd{k}.xVar;
                ylist = cpd{k}.data;
                errlist = cpd{k}.err;               
                errorbar(xlist+xShift,ylist(:,i),errlist(:,i),'-o','Color', hex2RGB(colorCodes{i}{k}), 'MarkerSize', 8, 'LineWidth', 1.5);
                ymin = [ymin min(ylist(:,i))];
                ymax = [ymax max(ylist(:,i))]; 
                xVarPlot_yMin = min(0,1.2.*min(ymin));
                xVarPlot_yMax = max(0.1,1.2.*max(ymax));
                xlist1 = cpd{1}.xVar;
                rep = cpd{k}.rep;
%                 text(min(xlist1) + 0.5.*range(xlist1), 0.4*(xVarPlot_yMax - xVarPlot_yMin)+0.05*k*(xVarPlot_yMax - xVarPlot_yMin), b.files{k}, 'FontSize', 13, 'HorizontalAlignment', 'right', 'Color',  hex2RGB(colorCodes{i}{k}), 'Interpreter', 'none')
                text(min(xlist1) + 0.25.*range(xlist1), 0.25*(xVarPlot_yMax - xVarPlot_yMin)+0.08*k*(xVarPlot_yMax - xVarPlot_yMin), [b.comments{k},' ',num2str(rep),' ','cycles'], 'FontSize', 13, 'HorizontalAlignment', 'right', 'Color',  hex2RGB(colorCodes{i}{k}), 'Interpreter', 'none')
            end
            xVarPlot_yMin = min(0,1.2.*min(ymin));
            xVarPlot_yMax = max(0.1,1.2.*max(ymax));
%             if i==3
%                  xVarPlot_yMax = 1;
%             end
            ylim([xVarPlot_yMin,xVarPlot_yMax]);
           
            if b.rewriteXvar
                xlabel('690 laser frequency (MHz)', 'FontSize', 14)
            else
%                 xlabel('Valon frequency (MHz)', 'FontSize', 14)
%                 xlabel('Interaction time (us)', 'FontSize', 14)
%                 xlabel('IR laser frequency 296388 (GHz) + x [MHz]', 'FontSize', 14)
                xlabel(b.xlabel, 'FontSize', 14)
            end
            ylabel('Ion counts per cycle', 'FontSize', 14)
            
            text(min(xlist) + 1.1.*range(xlist), 0.50*(xVarPlot_yMax - xVarPlot_yMin), [b.speciesNames{i},'^+'], 'FontSize', 16, 'HorizontalAlignment', 'right', 'Color',  hex2RGB(colorCodes{i}{1}))
            %             text(min(xlist) + 1.1.*range(xlist), 0.30*(xVarPlot_yMax - xVarPlot_yMin), [num2str(totCycleNumEff),' cycles'], 'FontSize', 16, 'HorizontalAlignment', 'right')
            box on
            grid on
            grid minor
            set(gca,'FontSize',16)
        end
    end      
end

function plot_background(cpd,b)
    h7 = figure(); 
    scrsz = get(groot, 'Screensize');
    set(h7, 'Position', [scrsz(3) + scrsz(3)/2 1 scrsz(3)/2 scrsz(4)])
    set(h7,'Color','w')
    species = b.species;
    pinks = {'#860454','#c51b7c','#dc75ab','#f0b7da','#FF00FF','#ff81c0'};
    greens = {'#276418','#4d921e','#7fbc42','#b9e084','#15b01a','#77AC30'};
    purples = {'#410149','#762a84','#9b6fac','#c1a5cd','#800080','#7e1e9c'};
    browns = {'#523107','#7f3c0a','#bf812c','#e2c17e','#a0522d','#A2142F'};
    pools = {'#003d2e','#01665e','#379692','#81cdc1','#13eac9','#04d8b2'};
    reds = {'#68011d','#b5172f','#d75f4e','#f7a580','#ff0000','#e50000'};
    blues = {'#062e61','#2265ad','#4295c1','#93c5dc','#40e0d0','#00FFFF'};
    b.species = [0, 0, 1, 1, 0, 1 , 1, 0];
    colorCodes = {reds,pools,purples,greens,pinks,blues,browns,pinks};
    for i = 1:length(b.speciesNames)
        if species(i)
            figure(h7)
            hold on
            if length(nonzeros(species)) <= 4
                subplot(length(nonzeros(species)),1,find(find(species == 1)==i))
                hold on
            else
                subplot(ceil(length(nonzeros(species))/2),2, find(find(species == 1)==i))
                hold on
            end
            
            ymin = [];
            ymax = [];
            for k = 1:length(cpd)
                xShift = b.xShift{k};
                xlist = cpd{k}{1}.xVar;
                ylist = cpd{k}{1}.data;
                errlist = cpd{k}{1}.err;        
                
                if k==2
                    ylist = cpd{k}{1}.data;
                    errlist = cpd{k}{1}.err;
                end
                errorbar(xlist+xShift,ylist(:,i),errlist(:,i),'-o','Color', hex2RGB(colorCodes{i}{2*k-1}), 'MarkerSize', 8, 'LineWidth', 1.5);
                ymin = [ymin min(ylist(:,i))];
                ymax = [ymax max(ylist(:,i))]; 
                xVarPlot_yMin = min(0,1.2.*min(ymin));
                xVarPlot_yMax = max(0.1,1.2.*max(ymax));
                xlist1 = cpd{k}{1}.xVar;
                rep = cpd{k}{1}.rep;
%                 text(min(xlist1) + 0.5.*range(xlist1), 0.4*(xVarPlot_yMax - xVarPlot_yMin)+0.05*k*(xVarPlot_yMax - xVarPlot_yMin), b.files{k}, 'FontSize', 13, 'HorizontalAlignment', 'right', 'Color',  hex2RGB(colorCodes{i}{k}), 'Interpreter', 'none')
                text(min(xlist1) + 0.25.*range(xlist1), 0.25*(xVarPlot_yMax - xVarPlot_yMin)+0.08*k*(xVarPlot_yMax - xVarPlot_yMin), [b.comments{k},' ',num2str(rep),' ','cycles'], 'FontSize', 13, 'HorizontalAlignment', 'right', 'Color',  hex2RGB(colorCodes{i}{2*k-1}), 'Interpreter', 'none')
            end
            xVarPlot_yMin = min(0,1.2.*min(ymin));
            xVarPlot_yMax = max(0.1,1.2.*max(ymax));
%             if i==3
%                  xVarPlot_yMax = 2;
%             end
            ylim([xVarPlot_yMin,xVarPlot_yMax]);
           
            if b.rewriteXvar
                xlabel('690 laser frequency 434920 GHz + x (MHz)', 'FontSize', 14)
            else
%                 xlabel('Valon frequency (MHz)', 'FontSize', 14)
%                 xlabel('Interaction time (us)', 'FontSize', 14)
%                 xlabel('IR laser frequency 296388 (GHz) +x [MHz]', 'FontSize', 14)
                 xlabel(b.xlabel, 'FontSize', 14)
            end
            ylabel('Ion counts per cycle', 'FontSize', 14)
            
            text(min(xlist) + 1.1.*range(xlist), 0.50*(xVarPlot_yMax - xVarPlot_yMin), [b.speciesNames{i},'^+'], 'FontSize', 16, 'HorizontalAlignment', 'right', 'Color',  hex2RGB(colorCodes{i}{1}))
            %             text(min(xlist) + 1.1.*range(xlist), 0.30*(xVarPlot_yMax - xVarPlot_yMin), [num2str(totCycleNumEff),' cycles'], 'FontSize', 16, 'HorizontalAlignment', 'right')
            box on
            grid on
            grid minor
            set(gca,'FontSize',16)
        end
    end  
    
    if 1 % for normalization
        h=figure(171)
        xref = cpd{1}{1}.xVar;
        yref = cpd{1}{1}.data(:,3);
        errref = cpd{1}{1}.err(:,3);
        
        
        ydiff = cpd{1}{1}.data(:,3)-cpd{2}{1}.data(:,3);
        err_diff = sqrt(cpd{1}{1}.err(:,3).^2+cpd{2}{1}.err(:,3).^2);
        ynorm = ydiff/max(ydiff);%./yref/max(ydiff./yref);
        enorm = err_diff/max(ydiff);
%         errorbar(xref,ynorm,enorm,'-^','markersize',10,'linewidth',2);
%         errorbar(sin(xref/37.1*pi/2).^2, ydiff,err_diff,'-^','markersize',10,'linewidth',2);
        errorbar(xref, ydiff,err_diff,'-o','markersize',10,'linewidth',2);
        hold on 
        ynorm1 = ydiff./yref/(max(ydiff./yref));
        enorm1 = ynorm1.*sqrt(err_diff.^2./(ydiff.^2)+errref.^2./(yref.^2))/(max(ydiff./yref));
%         plot(xref,ynorm1,'-o','markersize',10,'linewidth',2)
%         errorbar(xref,ynorm1,enorm1,'-o','markersize',10,'linewidth',2)
      
%         xlabel('Interaction time (us)')
%         xlabel('IR laser frequency 296388 (GHz) +x [MHz]', 'FontSize', 14)
        xlabel(b.xlabel, 'FontSize', 14)
%         xlabel('sin(t/t_{\pi}*\pi/2)^2', 'FontSize', 14)
        ylabel('Signal')
        box on; grid on; grid minor; 
        legend('y1-y2 ','(y1-y2)/yRb, normalized')
        set(gca,'Fontsize',16)
        
%         h=figure(172);
%         ratio = [0    0.1067    0.2456    0.4213    0.6462    0.8865    1.0000];
%         err_ratio = [0    0.0175    0.0199    0.0309    0.0355    0.0663    0.0587];
%         errorbar(ratio,ydiff,err_diff,err_diff,err_ratio,err_ratio,'-o','markersize',10,'linewidth',2);
%         xlabel('Normalized atom number');
%         ylabel('KRb ion counts')
%         box on; grid on; grid minor; 
%         set(gca,'Fontsize',16)
    end
end


function plot_overlap_species(cpd,b,ind)
    figure(2)
    pinks = {'#860454','#c51b7c','#dc75ab','#f0b7da','#FF00FF','#ff81c0'};
    greens = {'#276418','#4d921e','#7fbc42','#b9e084','#15b01a','#77AC30'};
    purples = {'#410149','#762a84','#9b6fac','#c1a5cd','#800080','#7e1e9c'};
    browns = {'#523107','#7f3c0a','#bf812c','#e2c17e','#a0522d','#A2142F'};
    pools = {'#003d2e','#01665e','#379692','#81cdc1','#13eac9','#04d8b2'};
    reds = {'#68011d','#b5172f','#d75f4e','#f7a580','#ff0000','#e50000'};
    blues = {'#062e61','#2265ad','#4295c1','#93c5dc','#40e0d0','#00FFFF'};
    b.species = [0, 0, 1, 1, 0, 1 , 1, 0];
    colorCodes = {reds,pools,purples,greens,pinks,blues,browns,pinks};
    for i = 1:length(b.speciesNames)
        if b.species(i)         
            xlist = cpd{ind}.xVar;
            ylist = cpd{ind}.data;
            errlist = cpd{ind}.err;               
            errorbar(xlist,ylist(:,i),errlist(:,i),'-o','Color', hex2RGB(colorCodes{i}{ind}), 'MarkerSize', 8, 'LineWidth', 1.5);
            hold on
            ymin =  min(ylist(:,i));
            ymax = max(ylist(:,i)); 
            xVarPlot_yMin = min(0,1.2.*min(ymin));
            xVarPlot_yMax = max(0.1,1.2.*max(ymax));
            xlist1 = cpd{1}.xVar;
%                 text(min(xlist1) + 0.5.*range(xlist1), 0.4*(xVarPlot_yMax - xVarPlot_yMin)+0.05*k*(xVarPlot_yMax - xVarPlot_yMin), b.files{k}, 'FontSize', 13, 'HorizontalAlignment', 'right', 'Color',  hex2RGB(colorCodes{i}{k}), 'Interpreter', 'none')
       

%             if i==3
%                  xVarPlot_yMax = 4;
%             end
            ylim([xVarPlot_yMin,xVarPlot_yMax]);
           
            if b.rewriteXvar
                xlabel('IR laser frequency (MHz)', 'FontSize', 14)
                xlabel(b.xlabel, 'FontSize', 14)
            else
                xlabel('Valon frequency (MHz)', 'FontSize', 14)
                xlabel(b.xlabel, 'FontSize', 14)
            end
            ylabel('Ion counts per cycle', 'FontSize', 14)
            
            text(min(xlist) + 1.1.*range(xlist), 0.50*(xVarPlot_yMax - xVarPlot_yMin), [b.speciesNames{i},'^+'], 'FontSize', 16, 'HorizontalAlignment', 'right', 'Color',  hex2RGB(colorCodes{i}{1}))
            %             text(min(xlist) + 1.1.*range(xlist), 0.30*(xVarPlot_yMax - xVarPlot_yMin), [num2str(totCycleNumEff),' cycles'], 'FontSize', 16, 'HorizontalAlignment', 'right')
            box on
            grid on
            grid minor
            set(gca,'FontSize',16)
        end
    end   
end

function color = hex2RGB(hex)
    color = sscanf(hex(2:end),'%2x%2x%2x',[1 3])/255;
end